PercolationAndCaprise Subroutine

SUBROUTINE PercolationAndCaprise  & 
!
(id, i, j, rain, vadose, pet, soilDepthCell, percolationcellRZ, &
    percolationcellTZ, runoffcell)

IMPLICIT NONE

! Arguments with intent(in):
INTEGER (KIND = short), INTENT(IN) :: id !!soil balance id
INTEGER (KIND = short), INTENT(IN) :: i,j
REAL (KIND = float), INTENT(IN) :: rain
TYPE (grid_real), INTENT(IN) :: vadose !!vadose zone depth
TYPE (grid_real), INTENT(IN) :: pet !!potential evapotranspiration

!Arguments with intent(out):
REAL (KIND = double), INTENT(OUT) :: soilDepthCell
REAL (KIND = float), INTENT(OUT) :: percolationcellRZ
REAL (KIND = float), INTENT(OUT) :: percolationcellTZ
REAL (KIND = double), INTENT(OUT) :: runoffcell

!Local declaration:
REAL (KIND = float) :: smAdjustedRZ !! used to prevent soil moisture going to zero
REAL (KIND = float) :: smAdjustedTZ !! used to prevent soil moisture going to zero
REAL (KIND = float) :: conductivityCellRZ !!partial saturation soil hydraulic 
                                        !conductivity of current cell [m/s]
REAL (KIND = float) :: conductivityCellTZ !!partial saturation soil hydraulic 
                                        !conductivity of current cell [m/s]
                                        !conductivity of current cell [m/s]
REAL (KIND = float) :: psiCell !!matric potential of current cell [m]
REAL (KIND = float) :: meanHydCond !!mean hydraulic conductivity used to 
                                   !compute capillary rise [m/s]
REAL (KIND = float) :: dsdt
!------------end of declaration------------------------------------------------

!calculate unsaturated hydraulic conductivity [m/s]

!root zone
IF (soilMoistureRZ % mat(i,j) <= thetar % mat(i,j)) THEN
    smAdjustedRZ = thetar % mat(i,j) + 0.001
ELSE
    smAdjustedRZ = soilMoistureRZ % mat(i,j) 
END IF

conductivityCellRZ = UnsHydCond (ksat = ksat % mat(i,j), &
                    theta = smAdjustedRZ, &
                    thetas = thetas % mat(i,j), &
                    thetar = thetar % mat(i,j), &
                    psdi = psdi % mat(i,j) ) 

!transmission zone
IF (soilMoistureTZ % mat(i,j) <= thetar % mat(i,j)) THEN
    smAdjustedTZ = thetar % mat(i,j) + 0.001
ELSE
    smAdjustedTZ = soilMoistureTZ % mat(i,j) 
END IF

conductivityCellTZ = UnsHydCond (ksat = ksat % mat(i,j), &
                    theta = smAdjustedTZ, &
                    thetas = thetas % mat(i,j), &
                    thetar = thetar % mat(i,j), &
                    psdi = psdi % mat(i,j) ) 

IF (id == LANDPLAIN ) THEN

         
    !set local soil depth
    !soilDepthCell = MIN (soildepth % mat(i,j), vadose % mat(i,j) )
    soilDepthCell = soildepth % mat(i,j)
    
    !compute root zone soil suction (m)
    psiCell = Psi (psic = psic % mat(i,j), theta = smAdjustedRZ, &
                    thetas = thetas % mat(i,j), &
                    thetar = thetar % mat(i,j), &
                    psdi = psdi % mat(i,j) ) 
				   
    !Compute harmonic mean among saturated conductivity 
    !at the interface with groundwater table and unsaturated 
    !conductivity of vadose zone
    meanHydCond = 2. * conductivityCellRZ * ksat % mat(i,j) / &
                  (conductivityCellRZ + Ksat % mat(i,j))
        

    !compute capillary rise and percolation 
    IF  ( vadose % mat (i,j) < soilDepthRZ % mat (i,j) )  THEN
        ! water table depth lies within the root zone or above ground
        saturatedByGroundwater = .TRUE.
    ELSE
        capRise % mat (i,j) =  meanHydCond * psiCell / vadose % mat (i,j)
        dsdt = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthRZ % mat(i,j) / dtSoilBalance
        percolationcellRZ = MIN ( conductivityCellRZ, dsdt)
            
        dsdt = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthTZ % mat(i,j) / dtSoilBalance
        percolationcellTZ = MIN ( conductivityCellTZ, dsdt) * &
                                percolationFactor % mat (i,j)
        
        saturatedByGroundwater = .FALSE.
    END IF
       
				    
ELSE !slope or channel cell
			soilDepthCell = soilDepth % mat(i,j)
			capRise % mat(i,j) = 0.
            dsdt = ( soilMoistureRZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthRZ % mat(i,j) / dtSoilBalance
            percolationcellRZ = MIN ( conductivityCellRZ, dsdt)
            
            dsdt = ( soilMoistureTZ % mat(i,j) - thetar % mat(i,j) ) * &
                      soilDepthTZ % mat(i,j) / dtSoilBalance
            percolationcellTZ = MIN ( ksat % mat(i,j), dsdt) * &
                                percolationFactor % mat (i,j)
END IF

 
RETURN
END SUBROUTINE PercolationAndCaprise